global nx ny nz dif wi w vex vey ky phi T T0 T0i Ti dyT...
    dy1p5 dx1p5 dyt1p5 dxt1p5 dzt1p5 dtx dty Ax A3 phi2 xx yy...
     vii dx dy dz dt ak2 nx0 ny0 nz0 phii pert

nx=nx0+2;
ny=ny0+2;
nz=nz0+2;
dx=alx/(nx0+1);
dy=aly/ny0;
dz=alz/nz0;
dxt1p5=tau/(2.*dx);
dyt1p5=tau/(2.*dy);
dzt1p5=tau/(dz*dz);

dx1p5=1./(2.*dx);
dy1p5=1./(2.*dy);
dtx=tau/(dx^2);
dty=tau/(dy^2);
dtdx=tau/dx;

%定义了一个对称的三对角矩阵 Ax,
%主对角线上的每个元素都是 2,次对角线（上边和下边紧邻主对角线的两条对角线）上的每个元素都是 -1
% Ax=2*eye(nx0)-diag(ones(nx0-1,1),1)-diag(ones(nx0-1,1),-1);
% %除以dx^2变成差分矩阵算子，再稀疏化
% Ax=Ax/dx^2;
% Ax=sparse(Ax);

% ny02=ny0/2;
% pi2=2.*pi;

% pi2y=pi2/aly;   %2pi/Ly
% vkysq=zeros(ny,1);
% vky=zeros(ny,1);
% for j=1:ny0/2+1
%     vky(j)=pi2y*(j-1);  %2pi/Ly*(j-1)
%     vkysq(j)=vky(j)^2;%(2pi/Ly*(j-1))^2
% end
% for j=ny0/2+2:ny0
%     jj=j-ny0;%负指标
%     vky(j)=pi2y*(jj-1);
%     vkysq(j)=vky(j)^2;
% end

% %vkyz = vky;
% vkysqz = vkysq;


% ak2=zeros(ny,1);
% ak2(1:ny0,1)=vkysqz(1:ny0);


% x2=zeros(nx,1);
% for i=1:nx
%     x2(i)=dx*(i-1)-alx/2;
% end
% %x2y = repmat(x2,[1,ny]);
% %vkyx = repmat(vky',[nx,1]);

% A4=kron(ak2(1:ny0,1),eye(nx0))+kron(ones(ny0,1),Ax);
% A4=sparse(A4);
% A3=sparse(nx0*ny0,nx0*ny0);
% for i=1:ny0
%     v1=sparse(i,i,1,ny0,ny0);
%     E=kron(v1,eye(nx0));
%     u=kron([zeros(1,i-1),1,zeros(1,ny0-i)],eye(nx0));
%     u=sparse(u);
%     incr=E*A4*u;
%     incr=sparse(incr);
%     A3=A3+incr;
% end
% save('A3.mat','A3');
% clear A4;
A3 = generate_poisson_matrix(nx0, ny0, alx,  aly); 


%% Initialization
wi=zeros(nx,ny);
Ti=zeros(nx,ny);
T0i=zeros(nx,ny);
vex=zeros(nx,ny);
vey=zeros(nx,ny);
dyT=zeros(nx,ny);
vii=zeros(nx,ny);
phii=zeros(nx,ny);
phi2=zeros(nx,ny);

% T0_fun = @(x,y)sinh(ky-ky*x)/sinh(ky);
T0_fun = @(x,y)1-x;

xx=0:dx:alx;yy=0:dy:aly+dy;zz = 0:dz:alz+dz;
[x,y] = ndgrid(xx,yy);
% filter=reshape(exp(-(xx-alx/2).^2/6),[nx,1])*ones(1,ny);

filty = [zeros(nx,1),ones(nx,9),zeros(nx,ny-21),ones(nx,9)];
filtx = [zeros(ny-2,1),ones(ny-2,9),zeros(ny-2,nx-19),ones(ny-2,9)]';
T0i(:,:)=T0_fun(x/alx,y);

sa=rand(nx,ny-2)-0.5;sv=fft(sa,[],2);
sa=sv.*filty;
wi(:,2:end-1)=ifft(sa,[],2);
fftw = fft(wi(:,2:end-1),[],1).*filtx;
wi(:,2:end-1)=ifft(fftw,[],1);

sa1=rand(nx,ny-2)-0.5;sv1=fft(sa1,[],2);
sa1=sv1.*filty;
Ti(:,2:end-1)=ifft(sa1,[],2);
fftT = fft(Ti(:,2:end-1),[],1).*filtx;
Ti(:,2:end-1)=ifft(fftT,[],1);



% K = (1/64)*ones(8);
% wi = conv2(wi,K,'same');
% Ti = conv2(Ti,K,'same');
% for i = 1:nx
%     wi(i,2:end-1) = smooth(wi(i,2:end-1),5);
%     Ti(i,2:end-1) = smooth(Ti(i,2:end-1),5);
% end



wi=sbcxnn3(wi*pert*0.01);wi=sbcy(wi);
Ti=sbcxnn3(Ti*pert);Ti=sbcy(Ti);

spe1d_fix_p_mod
%     phi0=sbcy(phi0);
%     wi0=sbcy(wi0);
%     phii=sbcy(phii);
%     phi=phii+phi0;
if restart == 0  
    fid=['initial'];
    save([saldir,fid]);

    pcolor(yy,xx,squeeze(wi(:,:))); 
    colorbar; colormap jet; shading interp;
    title('$$\omega$$','Interpreter','latex')
    ylabel('x')
    xlabel('y')
    print(gcf,'-dpng',[saldir1,'initial_omega'])
    close
    
    
    pcolor(yy,xx,squeeze(Ti(:,:))); 
    colorbar; colormap jet; shading interp;
    title('T','Interpreter','latex')
    ylabel('x')
    xlabel('y')
    print(gcf,'-dpng',[saldir1,'initial_T'])
    close
    
    pcolor(yy,xx,squeeze(phii(:,:))); 
    colorbar; colormap jet; shading interp;
    title('T0','Interpreter','latex')
    ylabel('x')
    xlabel('y')
    print(gcf,'-dpng',[saldir1,'initial_phi'])
    close
    
end
%% solve equations

w=wi;T=Ti;T0=T0i;

if (restart == 1)
   fid_start=['dat',sprintf('%4.4d',nts_start)];
   fid_start=['./data/',fid_start,'.mat'];
   load(fid_start)
end
for nt=1:nts
    for ntt=1:ntp
        sfield
        getdyT
        f=0.5;
        fi=0.5;
        sT(f,fi);
        sw(f,fi);
%         sT0(f,fi);
        spe1d_fix_p_mod
        %         phi=phi0+phii;
        sfield
        getdyT
        f=1.0;
        fi=0.0;
        sT(f,fi);
        sw(f,fi);
%         sT0(f,fi);
        spe1d_fix_p_mod
        %         phi=phi0+phii;
        wi=sbcy(wi);wi=sbcxnn3(wi);
        Ti=sbcy(Ti);Ti=sbcxnn3(Ti);
%         T0i=sbcT0(T0i);
        done=isfinite(wi(2,2));
        if done==0
            error('nan')
        end
    end 
    if (restart==1)
        disp(nt+nts_start);
%         sprintf(['Time is t = ',num2str(nt+nts_start)])
        fid=['dat',sprintf('%4.4d',nt+nts_start)];
    else
    sprintf(['Time is t = ',num2str(nt)])
    fid=['dat',sprintf('%4.4d',nt)];
    end
    save([saldir,fid],'wi','phii','T0i','Ti','vex','vey')% 'phi2','tei')%,'pii', 'pei','vii')%,'laplace2_pi')
end
